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Abstract 



Euler-Lagrange variational principle is used to obtain analytical and numerical 
flow relations in cylindrical tubes. The method is based on minimizing the total 
stress in the flow duct using the fluid constitutive relation between stress and rate 
of strain. Newtonian and non-Newtonian fluid models; which include power law, 
Bingham, Herschel-Bulkley, Carreau and Cross; are used for demonstration. 

Keywords: Euler-Lagrange variational principle; fluid mechanics; generalized 
Newtonian fluid; capillary flow; pressure-flow rate relation; Newtonian; power law; 
Bingham; Herschel-Bulkley; Carreau; Cross. 



1 INTRODUCTION 



6 



1 Introduction 

Several methods are in use to derive relations between pressure, p, and volumet- 
ric flow rate, Q, in tubes and conduits. These methods include the application 
of first principles of fluid mechanics with utilizing the fluid basic properties [1- 
3], the use of Navier-Stokes equations [4], the lubrication approximation [5], and 
Weissenberg-Rabinowitsch-Mooney relation [1-3, 6]. Numerical methods related to 
these analytical formulations, such as finite element and similar meshing techniques 
[7], are also in use when analytical expressions are not available. 

However, we are not aware of the use of the Euler-Lagrange variational principle 
to derive p-Q relations in general and in capillaries in particular despite the fact 
that this principle is more intuitive and natural to use as it is based on a more 
fundamental physical principle which is minimizing the total stress combined with 
the utilization of the fluid constitutive relation between stress and rate of strain. 

The objective of this paper is to outline this method demonstrating its appli- 
cation to Newtonian and some time-independent non-Newtonian fluids and fea- 
turing its applicability numerically as well as analytically. In the following, we 
assume a laminar, axi-symmetric, incompressible, steady, viscous, isothermal, fully- 
developed flow for generalized Newtonian fluids moving in cylindrical tubes where 
no-slip at wall condition [8] applies and where the flow velocity profile has a sta- 
tionary derivative point at the middle of the tube (r = 0) meaning the profile has 
a blunt rounded vertex. 

2 Method 

The constitutive relation for generalized Newtonian fluids in shear flow is given by 



t = m 



(1) 
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where r is the shear stress, 7 is the rate of shear strain, and y is the shear viscosity 
which generally is a function of the rate of shear strain. It is physically intuitive 
that the flow velocity profile in a tube (or in a flow path in general) will adjust 
itself to minimize the total stress which is given by 



. ' f R dr , f R d . , , r f da dj\ , 

Tt= I dT= L Tr dT= L Tr^ l)dr= L VTr + »Tr ]dr (2) 



where r t is the total stress, r c and t w are the shear stress at the tube center and tube 
wall respectively, and R is the tube radius. The total stress, as given by Equation 
2, can be minimized by applying the Euler-Lagrange variational principle which, 
in its most famous form, is given by 



df d fdf\ 

(3) 



dy dx \dy' 
where 



d\L dj df d ( d[L dj\ , 

x = r, y = y, / = 7— , and — = — 7— + \i— = fi 4 

dr dr ay oy \ dr dr 1 



However, to simplify the derivation we use here another form of the Euler-Lagrange 
principle which is given by 



±(f- >?l\_?f = Q (5) 
dx \ dy' J dx 

that is 



d ( dfi drt dy\ d ( dy, dry\ 

i~r + y~r - y~r } - ^ [i-r + y~r } =® (6) 



dr V dr dr dr J dr V dr dr 



i.e. 
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A (ryii) - — (r^E + fj,^l\ = (7) 

dr \ dr J dr \ dr dr J 

Since ordinary derivative is a special case of partial derivative, we can write this 
equation as 



d f dji d\i dj 



dr dr ^ dr ^ dr 



S ^')=0 (9) 



that is 



dr \ dr 

In the following sections the use of this equation will be demonstrated to derive 
p-Q flow relations for generalized Newtonian fluids. 

2.1 Newtonian 

For Newtonian fluids, the viscosity is constant that is 

/I = fJL (10) 

and hence Equation 9 becomes 



dr \ dr 

On integrating once with respect to r we obtain 



C> ^$)=o (ID 



dr 

where A is a constant. Hence 



= A (12) 
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7 = — (Ar + B) (13) 

where .£> is another constant. Now from the two boundary conditions at r = and 
r = R, A and B can be determined, that is 

7 (r = 0)=0 5 = (14) 

and 

, n , r w PR AR P , . 

7 ( r = R) = JUL = — — = A = —— 15 

where r w is the shear stress at the tube wall, P is the pressure drop across the tube 
and L is the tube length. Hence 

7 (r) = ^-r (16) 

On integrating this with respect to r the standard Hagen-Poiseuille parabolic 
velocity profile is obtained, that is 



v (r) 



/[ dv f f P P 

dv = —dr = / 7<ir = / — — rdr = — — r 2 + D (17) 
J dr J J 2L/i 4L/i 



where v(r) is the fluid axial velocity at r and D is another constant which can be 
determined from the no-slip at wall boundary condition, that is 



PR 2 

v ( r = R ) = ^ D = -^_ (18) 

4L/i 

that is 
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where the minus sign at the front arises from the fact that the pressure gradient is 
opposite in direction to the flow velocity vector. The volumetric flow rate will then 
follow by integrating the flow velocity profile with respect to the cross sectional 
area, that is 

r R nP r R tt pr 4 

Q= \v\ 2nrdr = ^f- / (R 2 - r 2 ) rdr = (20) 
Jo 2L/i D J 8Lfi Q 

which is the well-known Hagen-Poiseuille flow relation. 
2.2 Power Law 

For power law fluids, the viscosity is given by 



H = k^ 1 (21) 
On applying Euler-Lagrange variational principle (Equation 9) we obtain 



d ( kf f-^\ =0 (22) 



dr \ dr 
On integrating once with respect to r we obtain 



kl n-i<tL = A (23) 
dr 

On separating the two variables and integrating both sides we obtain 



7 =W£(i4r + fl) (24) 

where A and B are constants which can be determined from the two boundary 
conditions, that is 



7 (r = 0)=0 B = 



(25) 
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and 




7 ( r = R) = W = {/^ = p/JaR A = — (26) 

1 ' w 2LA; V fc 2nL v ; 



^ „l/n 



and therefore 



Tiimr (27) 

On integrating this with respect to r the power law velocity profile is obtained, 
that is 



v(r)= ( dv = I —dr = [^dr= [ \ —r^dr = \ — r 1+1 / n + D 

y ' J J dr J 1 J V 2kL n + 1 V 2kL 

(28) 

where D is another constant. From the no-slip at wall boundary condition 



«(' = *>=<> ^—i^TVW^" (29) 



that is 



"M = ^iV2^( fll+1/ "- rl+1/ ") (30) 

The volumetric flow rate will then follow by integrating the velocity profile with 
respect to the cross sectional area, that is 



Q= f R \v\ 2-nrdr = — [j^r f (R l+l / n - r l+l / n ) rdr 
io n + 1 V 2kL J K > 



"" P -i?3+l/n (31) 



3n + 1 V 2kL 

which is mathematically equivalent to the expressions derived in [1-3, 5] using other 
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methods. 

3 Numerical Implementation 

For generalized Newtonian fluids with complex constitutive relations, it may be 
very difficult, or even impossible, to obtain a flow analytical solution from the 
Euler-Lagrange principle. In this case, the variational method can be used as a 
basis for a numerical method by employing Equation 9 to obtain the rate of shear 
strain as an explicit or implicit function of r which is then numerically solved 
and integrated to obtain the flow velocity profile which, in its turn, is numerically 
integrated to obtain the p-Q relation. For the fluids which have an explicit relation 
between the rate of strain and radius, such as Newtonian and power law fluids 
(refer to Equations 16 and 27), the rate of strain can be computed directly for each 
r. However, for the fluids which have no such explicit relation, such as Bingham, 
Herschel-Bulkley, Carreau and Cross (refer to Equations 33, 36, 39 and 46), a 
numerical solver, based for instance on a bisection method, is required to find the 
rate of strain as a function of radius. A numerical integration scheme; such as 
midpoint, trapezium or Simpson rule; can then be utilized to integrate the strain 
rate with respect to radius to obtain the velocity profile in the first stage, and 
to integrate the velocity profile with respect to the cross sectional area to get the 
volumetric flow rate in the second stage. 

The constant of the first integration (strain rate with respect to radius) is in- 
corporated within a boundary condition by starting at the wall with zero velocity 
(v = 0); the velocity growth at the next inner ring of the tube cross section, 
obtained from numerically integrating the shear rate over radius, is then added 
incrementally to the velocity of the neighboring previous outer ring to obtain the 
velocity at the inner ring. The volumetric flow rate is then computed by multiply- 
ing the velocity of the ring with its cross sectional area and adding these partial 
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flow rate contributions to obtain the total flow rate. This method is applied, for 
the purpose of test and validation, to the Newtonian and power law fluids, for 
which analytical solutions are available, and the numerical results were compared 
to these analytical solutions. A sample of these comparisons are provided in Figures 
1 and 2. As seen, the numerical solutions converge fairly quickly to the analytical 
solutions; hence confirming the reliability of this numerical method and its theo- 
retical foundations. It should be remarked that the numerical results presented in 
the following sections were obtained by using three numerical integration schemes: 
midpoint, trapezium and Simpson rules. In all cases the three schemes converged 
to the same value although with different convergence rate. For the fluids with an 
implicit 7-r relation, a bisection numerical solver was used to obtain 7 as a function 
of r. 
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Figure 1: Q versus P plot for numeric solutions of a typical Newtonian fluid 
with fi Q = 0.005 Pa.s, flowing in a tube with L = 0.1 m and R = 0.01 m for 
r-discretization N r = 4, N r = 6 and N r > 100 alongside the analytical solution 
(Equation 20). 
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Figure 2: Q versus P plot for numeric solutions of a typical power law fluid with 
n = 0.9 and k = 0.005 Pa.s™, flowing in a tube with L = 0.1 m and R = 0.01 m 
for r-discretization N r = 4, N r = 6 and N r > 100 alongside the analytical solution 
(Equation 31). 

4 Yield Stress Fluids 

Experimenting with the use of Euler-Lagrange principle on different fluids, we 
tested this method on some yield stress fluids, specifically Bingham plastic and 
Herschel-Bulkley models, despite our awareness of the limitation of this method 
and its restriction to fluids which makes it inapplicable to yield stress materials 
since the solid-like plug flow at the center of the tube invalidates this assumption. 
For Bingham fluids, the viscosity is given by [1, 2, 9, 10] 



where r Q is the yield stress and C is the fluid consistency factor. On applying Euler- 
Lagrange variational principle (Equation 9) and following the method outlined in 
the Newtonian and power law fluid sections we obtain 



(32) 



7 
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T \nj + Cj = Ar + B (33) 

where A and B are the constants of integration. As seen in the last equation, the 
boundary condition at r = cannot be used to find B because 7 = is a singularity 
point. We therefore followed the non-yield stress fluid style and arbitrarily set 
B = 0. On applying the other boundary condition at r — R, we obtain 



For Herschel-Bulkley fluids, the viscosity is given by [1, 2, 9, 10] 



^ = _£ + C^- 1 (35) 
7 

Following a similar approach to that outlined in the Bingham part, the following 
relation was obtained 



C 

r hi7+ -7™ = Ar (36) 
n 

where 



These 7-r relations (i.e. Equation 33 for Bingham and Equation 36 for Herschel- 
Bulkley) were then solved for 7 at each r and numerically integrated to obtain the 
flow velocity profile first and volumetric flow rate second. The numerical results 
were interesting as the low yield stress fluids converged correctly to the analytic 
solution (refer to Figures 3 and 4) especially at high flow rates, while the high 
yield stress fluids diverged with finer discretization (refer to Figures 5 and 6). This 
can be explained by the occurrence of plug flow at the middle of the tube which 
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is considerable in the case of high yield stress fluids. The minor deviation of the 
low yield stress fluids at low pressures highlights this fact since the plug flow effect 
diminishes as the pressure and flow rate increase. The failure of this approach for 
high yield stress fluids is simply because our model applies to fluids; and yield stress 
materials prior to yield at plug area behave like solids. The convergence of the low 
yield stress fluids is due to the fact that negligible plug occurs at the middle of the 
tube especially at high pressures. 
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10000 



Figure 3: Q versus P plot for numeric solutions of a typical Bingham fluid above 
the yield point with C = 0.01 Pa.s and r Q = 1.0 Pa flowing in a tube with L = 0.1 m 
and R = 0.01 m for r-discretization N r = 4, N r = 14 and N r > 100 alongside the 
analytic solution [1, 2]. 



5 Carreau and Cross Fluids 

The Euler-Lagrange variational method was also applied to Carreau and Cross 
fluids for which no analytical expressions are available; the details are outlined in 
the following. 

For Carreau fluids, the viscosity is given by [1, 11-14] 
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Figure 4: Q versus P plot for numeric solutions of a typical Herschel-Bulkley fluid 
above the yield point with C = 0.01 Pa.s™, n = 0.8 and r D = 1.0 Pa flowing in a 
tube with L = 0.1 m and R = 0.01 m for r-discretization N r = A, N r = 10 and 
N r > 100 alongside the analytic solution [1, 2]. 
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Figure 5: Q versus P plot for numeric solutions of a typical Bingham fluid above the 
yield point with C = 0.01 Pa.s and r G = 10.0 Pa flowing in a tube with L = 0.1 m 
and R = 0.01 m for r-discretization N r — 4, N r — 10 and iV r > 100 alongside the 
analytic solution [1, 2]. 
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Figure 6: Q versus P plot for numeric solutions of a typical Herschel-Bulkley fluid 
above the yield point with C = 0.01 Pa.s n , n = 0.8 and r D = 10.0 Pa flowing in a 
tube with L = 0.1 m and R = 0.01 m for r-discretization N r = 4, N r = 10 and 
N r > 100 alongside the analytic solution [1, 2]. 



A* = Hoc + {Ho — A*oo) [l + (A7) 2 ] n (38) 

where [i a is the zero-shear viscosity, [i^ is the infinite-shear viscosity, A is a time 
constant, and n is the flow behavior index. On applying Euler-Lagrange variational 
principle (Equation 9) and following the derivation, as outlined in § 2.1 and 2.2, 
we obtain 

/l 1 — n 3 \ 

/ioo7 + (^o-/ioo)7 2Fi f 2' — ^ — ' 2 ; _A2 ^ 2 ) = Ar + B ( 39 ) 

where 2 -Pi is the hypergeometric function, and A and S are the constants of inte- 
gration. Now from the two boundary conditions at r = and r = R, A and B can 
be determined, that is 



7 (r = 0)=0 B = 



(40) 
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and 



/l 1 — n 3 \ 

7 (?" = #) = 7™ ^oclw + {^o-^oo)l W 2Fi[-,^—\-\-\ 2 -i 2 w \=AR 

(41) 

where 7^ is the shear rate at the tube wall. Now, by definition we have 

t w = Hwlw (42) 

that is 

^ = /ioo + {fl>o - fJ-oo) [l + (>^lw) 2 ] {n 1)12 Iw (43) 

From the last equation, j w can be obtained numerically by a numerical solver, 
based for example on a bisection method, and hence from Equation 41 A is obtained 

a Hoolw ~\~ (A*o A* 00) 7w 2-^1 (2) 2 ' 2 ' ^ Tw) / « < \ 

Several types of Carreau fluids were used for testing the model and its numerical 
implementation; a sample of these tests is presented in Figure 7 where the flow 
did converge for a radius discretization N r > 50. A sample of the flow velocity 
profile across the tube for the fluid of Figure 7 at a typical flow condition is also 
presented in Figure 8. These figures are qualitatively correct despite the fact that 
no analytical solution is available to fully validate the results. 

For Cross fluids, the viscosity is given by [15] 

. A*0 A* OO / , r\ 

^ = ^ + 1 + (A7) ™ ( 45 ) 

where \x Q is the zero-shear viscosity, floo is the infinite-shear viscosity, A is a time 
constant, and m is an indicial parameter. Following a similar derivation method 
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Figure 7: Q versus P plot for numeric solutions of a typical Carreau fluid with 
n = 0.75, fi = 0.05 Pa.s, yUoo = 0.001 Pa.s, and A = 1.0 s flowing in a tube with 
L = 0.5 m and R = 0.05 m for r-discretization N r = 10 and N r > 50. The numeric 
solutions were obtained using the real part of the hypergeometric function 2 Fi in 
Equations 39 and 44. 
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Figure 8: Flow velocity profile for the Carreau fluid of Figure 7 at a typical flow 
condition where the radius is scaled to unity and the speed is scaled to its maximum 
value at the tube center. 



to that outlined in Carreau, we obtain 
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/ioo7 + (fJ>o ~ Hoo) 7 2-Pl 1 



1 m + 1 



, -A 7 



Mi 



(46) 



77 i 



where 



.4 



/ioo7™ + (A*o - A*oo) 7w 2-^1 (l, 



m 1 m. 



1_. m+1 . 




(47) 



with 7^ being obtained numerically as outlined in Carreau. 

Several types of Cross fluids were used for testing the model and its numerical 
implementation; a sample of these tests is presented in Figure 9 where the flow 
did converge for a radius discretization A^ r > 50. A sample of the flow velocity 
profile across the tube for the fluid of Figure 9 at a typical flow condition is also 
presented in Figure 10. These figures are qualitatively correct despite the fact that 
no analytical solution is available to fully validate the results. The features in 
these results are similar to those observed in Carreau due to the strong similarities 
between these two fluids. 
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Figure 9: Q versus P plot for numeric solutions of a typical Cross fluid with 
n = 0.25, fi = 0.05 Pa.s, yUoo = 0.001 Pa.s, and A = 1.0 s flowing in a tube with 
L = 0.5 m and R = 0.05 m for r-discretization N r = 10 and N r > 50. The numeric 
solutions were obtained using the real part of the hypergeometric function 2 Fi in 
Equations 46 and 47. 
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Figure 10: Flow velocity profile for the Cross fluid of Figure 9 at a typical flow 
condition where the radius is scaled to unity and the speed is scaled to its maximum 
value at the tube center. 
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6 Conclusions 

In this article we outlined a method based on Euler-Lagrange variational principle 
which minimizes the total stress to obtain analytical and numerical flow relations for 
generalized Newtonian fluids in tubes and conduits in general. The method can be 
used in conjunction with numerical integration to obtain numerical solutions when 
analytical integration of the basic equations derived from the variational principle 
is difficult or impossible to obtain. The method was validated analytically and 
numerically for Newtonian fluids as well as a number of time-independent non- 
Newtonian fluids. 

The main advantages of this method are simplicity ease of implementation and 
rapid convergence to a solution which, for all practical purposes, is identical to 
the analytical solution. This convergence can be easily verified from two or more 
successive r-discretization schemes being converged to the same solution. 

The method can be used to obtain flow relations for complex fluids for which 
no analytical flow expressions have been derived from other methods due to math- 
ematical difficulties, such as Carreau, Carreau-Yasuda and Cross. This method, 
when implemented numerically in the case of analytical difficulties, is more accurate 
and suitable than the use of empirical relations or numerical meshing techniques. 

Numerical experiments were performed on Bingham and Herschel-Bulkley yield 
stress fluids to test the robustness of this method which is based on the assumption 
of fluidity. Interestingly, the method converged correctly to the analytical solution 
for low yield stress fluids although it did diverge for high yield stress fluids. The 
obvious reason which can explain these observations is that negligible plug flow 
occurs at the middle of the tube in the first case especially at high flow rates, 
while considerable plug flow occurs in the second case which invalidates the basic 
assumption of fluidity that this model relies upon. 

A preliminary investigation of the applicability of this method to Carreau and 
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Cross fluids has been conducted and presented in this article. More serious in- 
vestigations related to these fluids and other complex fluids by employing this 
variational technique are planned for the future. 
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Nomenclature 

7 shear rate (s _1 ) 

7 TO shear rate at tube wall (s _1 ) 

A time constant (s) 

fi fluid dynamic shear viscosity (Pa.s) 

\i zero-shear viscosity (Pa.s) 

/ioo infinite-shear viscosity (Pa.s) 

t shear stress (Pa) 

r c shear stress at tube center (Pa) 

t yield stress (Pa) 

r t total shear stress (Pa) 

r w shear stress at tube wall (Pa) 

C consistency factor in Bingham and Herschel-Bulkley models (Pa.s n ) 

2-Fi hypergeometric function 

k consistency factor in power law model (Pa.s™) 

L tube length (m) 

m indicial parameter in Cross model 

n flow behavior index 

iV r number of elements in radius discretization 

p pressure (Pa) 

P pressure drop (Pa) 

Q volumetric flow rate (m 3 .s _1 ) 

r radius (m) 
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R tube radius (m) 

v axial fluid velocity (m.s -1 ) 
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